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Abstract 

A  comprehensive  non-isothermal,  three-dimensional  computational  model  of  a  polymer  electrolyte  membrane  (PEM)  fuel  cell  has  been 
developed.  The  model  incorporates  a  complete  cell  with  both  the  membrane-electrode-assembly  (MEA)  and  the  gas  distribution  flow 
channels.  With  the  exception  of  phase  change,  the  model  accounts  for  all  major  transport  phenomena. 

The  model  is  implemented  into  a  computational  fluid  dynamics  code,  and  simulations  are  presented  with  an  emphasis  on  the  physical  insight 
and  fundamental  understanding  afforded  by  the  detailed  three-dimensional  distributions  of  reactant  concentrations,  current  densities, 
temperature  and  water  fluxes.  The  results  show  that  significant  temperature  gradients  exist  within  the  cell,  with  temperature  differences  of 
several  degrees  K  within  the  MEA.  The  three-dimensional  nature  of  the  transport  is  particularly  pronounced  under  the  collector  plates  land  area 
and  has  a  major  impact  on  the  current  distribution  and  predicted  limiting  current  density.  ©  2002  Elsevier  Science  B.V.  All  rights  reserved. 
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1.  Introduction 

Fuel  cells  convert  the  chemical  energy  of  hydrogen  and 
oxygen  directly  into  electricity.  Their  high  efficiency  and 
low  emissions  have  made  them  a  prime  candidate  for 
powering  the  next  generation  of  electric  vehicles,  and  their 
modular  design  and  the  prospects  of  micro-scaling  them 
have  gained  the  attention  of  cellular  phone  and  laptop 
manufacturers.  Their  scalability  and  relative  flexibility  in 
terms  of  fuel  makes  them  prime  candidates  for  a  variety  of 
stationary  applications  including  distributed  residential 
power  generation.  The  basic  structure  and  operation  prin¬ 
ciple  of  the  polymer  electrolyte  membrane  (PEM)  fuel  cell 
considered  here  are  illustrated  in  Fig.  1. 

The  polymer  electrolyte  consists  of  a  perfluorinated  poly¬ 
mer  backbone  with  sulfonic  acid  side  chains.  When  fully 
humidified,  this  material  becomes  an  excellent  protonic 
conductor.  The  membrane,  the  catalyst  (platinum  supported 
on  carbon  particle)  and  the  two  electrodes  (teflonated  porous 
carbon  paper  or  cloth)  are  assembled  into  a  sandwich 
structure  to  form  a  membrane-electrode-assembly  (MEA). 
The  MEA  is  placed  between  two  graphite  bipolar  plates  with 
machined  groves  that  provide  flow  channels  for  distributing 
the  fuel  (hydrogen)  and  oxidant  (oxygen  from  air). 
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The  hydrogen  rich  fuel  fed  to  the  anode  side  diffuses 
through  the  porous  gas-diffusion  electrode  (GDE).  At  the 
catalyst  layer,  the  hydrogen  splits  into  hydrogen  protons  and 
electrons  according  to: 

2H2  =>  4H+  +  4e_  (1) 

Driven  by  an  electric  field,  the  H+  ions  migrate  through  the 
polymer  electrolyte  membrane.  The  oxygen  in  the  cathode 
gas  stream  diffuses  through  the  towards  the  catalyst  interface 
where  it  combines  with  the  hydrogen  protons  and  the 
electrons  to  form  water  according  to: 

02  +  4H+  +  4e"  =+  2H20  (2) 

The  overall  reaction  is  exothermic  and  can  be  written  as: 

2H2  +  02  =+  2H20  +  electricity  +  heat  (3) 

Several  coupled  fluid  flow,  heat  and  mass  transport  processes 
occur  in  a  fuel  cell  in  conjunction  with  the  electrochemical 
reaction.  These  processes  have  a  significant  impact  on  two 
important  operational  issues:  (i)  thermal  and  water  manage¬ 
ment;  (ii)  mass  transport  limitations.  Water  management 
ensures  that  the  PEM  remains  fully  hydrated  to  maintain 
good  ionic  conductivity  and  performance.  Water  content  of 
the  membrane  is  determined  by  the  balance  between  water 
production  and  three  water  transport  processes:  electro- 
osmotic  drag  of  water,  associated  with  proton  migration 
through  the  membrane  from  the  anode  to  the  cathode  side; 
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Fig.  1.  Schematic  of  a  PEM  fuel  cell. 


back  diffusion  from  the  cathode;  and  diffusion  of  water  to/ 
from  the  oxidant/fuel  gas  streams.  Without  control,  an 
imbalance  between  production  and  removal  rates  of  water 
can  occur.  This  results  in  either  dehydration  of  the  mem¬ 
brane,  or  flooding  of  the  electrodes;  both  phenomena  have  a 
very  detrimental  effect  on  performance  and  fuel  cells  have  to 
be  carefully  designed  to  avoid  the  occurrence  of  either 
phenomenon.  Thermal  management  is  required  to  remove 
the  heat  produced  by  the  electrochemical  reaction  (up  to 
~50%  of  the  energy  produced  during  high  power  density 
operation)  in  order  to  prevent  drying  out  of  the  membrane 
and  excessive  thermal  stresses  that  may  result  in  rupture  of 
the  membrane.  The  small  temperature  differentials  between 
the  fuel  cell  stack  and  the  operating  environment  make 
thermal  management  a  challenging  problem  in  PEMFCs. 

Because  of  the  highly  reactive  environment  of  a  fuel  cell  it 
is  not  possible  to  perform  detailed  in  situ  measurements 
during  operation.  Such  information  has  been  sought  through 
modelling  and  simulation  in  order  to  improve  understanding 
of  water  and  species  transport,  optimize  thermal  manage¬ 
ment  and  shorten  the  design  and  optimization  cycles.  Mod¬ 
elling  of  fuel  cells  is  challenging,  because  the  processes 
involve  multi-component,  multi-phase,  and  multi-dimen¬ 
sional  flow,  heat  and  mass  transfer  with  electrochemical 
reactions,  all  occurring  in  irregular  geometries  including 
porous  media.  Numerous  authors  have  developed  fuel  cell 
models  accounting  for  various  physical  processes.  The  most 
prominent  earlier  works  stem  from  Bernardi  and  Verbrugge 
[3,4]  and  Springer  et.  al.  [14],  who  developed  one-dimen¬ 
sional,  isothermal  models  of  the  MEA.  Fuller  and  Newman 
[8]  published  a  quasi  two-dimensional  model  based  on 
concentration  theory.  The  work  by  Nguyen  and  White 
[12],  and  Yi  and  Nguyen  [21]  was  two-dimensional  in 
nature,  but  the  GDEs  were  omitted,  assuming  “ultrathin” 
electrodes.  The  importance  of  accounting  for  temperature 
gradients  in  fuel  cells  modelling  was  demonstrated  in  the 
work  of  Woehr  et.  al.  [20]  and  Djilali  and  Lu  [7].  The 
important  issue  of  water  flooding  was  addressed  by  Baschuk 
and  Li  in  a  recent  one-dimensional  model  [2]. 

Earlier  models  were  primarily  analytic  and  required  a 
number  of  simplifications  due  to  the  limitations  of  the 
numerical  techniques.  More  recently,  a  general  trend  can 


be  observed  to  apply  the  methods  of  computational  fluid 
dynamics  to  fuel  cell  modelling.  Gurau  et  al.  [9]  published  a 
fully  two-dimensional  model  of  a  whole  fuel  cell,  i.e.  two 
gas-flow  channels  separated  by  the  MEA.  Um  et  al.  [18]  and 
Wang  et  al.  [19]  have  developed  a  similar  model  and 
included  two  phase  flow.  However,  the  underlying  assump¬ 
tion  was  isothermal  behaviour,  which  is  a  serious  modelling 
limitation  as  we  will  see  later.  The  local  temperature  dis¬ 
tribution  has  a  significant  impact  on  the  amount  of  water  that 
undergoes  phase  change,  and  therefore  the  isothermal 
assumption  can  lead  to  results  that  are  not  physically 
representative  when  phase  change  is  accounted  for.  Finally 
as  a  result  of  the  architecture  of  a  cell,  the  transport 
phenomena  in  a  fuel  cell  are  inherently  three-dimensional, 
but  no  models  have  yet  been  published  to  address  this. 

In  this  paper  we  address  simultaneously  the  need  to 
account  for  thermal  gradients  and  multi-dimensional  trans¬ 
port  using  a  computational  fluid  dynamics  based  approach 
that  couples  convective  transport  in  the  gas-flow  channels 
with  transport  and  electrochemistry  in  the  MEA. 


2.  Model  description 

The  PEM  fuel  cell  model  presented  here  is  a  comprehen¬ 
sive  three-dimensional,  non-isothermal,  steady-state  model 
providing  a  detailed  description  of  the  following  transport 
phenomena: 

•  multi-component  flow; 

•  convective  heat  and  mass  transport  in  the  flow  channels; 

•  diffusion  of  reactants  through  porous  electrodes; 

•  electrochemical  reactions; 

•  migration  of  H+  protons  through  the  membrane; 

•  transport  of  water  through  the  membrane; 

•  transport  of  electrons  through  solid  matrix; 

•  conjugate  heat  transfer. 

The  equations  governing  these  processes  include  the  full 
mass  and  momentum  conservation  equations  (Navier- 
Stokes  equations)  governing  fluid  flow,  the  species  conser¬ 
vation  and  energy  equations  and  four  additional  phenom¬ 
enological  equations  tailored  to  account  for  processes 
specific  to  fuel  cells  [3]: 

•  the  Stefan-Maxwell  equations  for  multi-species  diffu¬ 
sion; 

•  the  Nernst-Planck  equation  for  the  transport  of  protons 
through  the  membrane; 

•  the  Butler- Volmer  equation  for  electrochemical  kinetics 
and; 

•  the  Schlogl  equation  for  the  transport  of  liquid  water 
through  the  membrane. 

These  equations  and  appropriate  boundary  conditions 
were  implemented  in  their  three-dimensional  form  into 
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the  commercial  computational  fluid  dynamics  code  CFX-4.3 
(AE A  Technology);  the  implementation  required  the  devel¬ 
opment  of  an  extensive  suite  of  user  subroutines.  Custo¬ 
mized  iterative  procedures  were  also  implemented  to  ensure 
effective  coupling  between  the  electrochemistry  and  the 
fluid  transport  processes. 

2.7.  Model  assumptions 

A  complete  fuel  cell  is  an  extremely  complex  system 
involving  both  microscale  and  macroscale  geometric  fea¬ 
tures  and  transport  processes.  In  order  to  devise  a  numeri¬ 
cally  tractable  three-dimensional  model  of  a  complete  cell  it 
is  necessary  to  invoke  a  number  of  simplifying  assumptions. 
The  most  important  ones  are: 

(1)  the  fuel  cell  operates  under  steady-state  conditions; 

(2)  all  gases  are  assumed  to  be  compressible  ideal  gases, 
fully  saturated  with  water  vapour; 

(3)  the  flow  in  the  channels  is  considered  laminar; 

(4)  the  membrane  is  assumed  to  be  fully  humidified  so 
that  the  electronic  conductivity  is  constant  and  no 
diffusive  terms  have  to  be  considered  for  the  liquid 
water  flux; 

(5)  since  it  was  determined  in  an  earlier  study  [4]  that 
cross-over  of  reactant  gases  can  be  neglected,  the 
membrane  is  currently  considered  impermeable  for 
the  gas  phase; 

(6)  the  water  in  the  pores  of  the  diffusion  layer  is 
considered  separate  from  the  gases  in  the  diffusion 
layers,  i.e.  no  interaction  between  the  gases  and  the 
liquid  water  exists; 

(7)  the  product  water  is  assumed  to  be  in  the  liquid  phase; 

(8)  Ohmic  heating  in  the  collector  plates  and  in  the  GDEs 
is  neglected  due  to  their  high  conductivity; 

(9)  heat  transfer  inside  the  membrane  is  accomplished  by 
conduction  only,  i.e.  the  enthalpy  carried  by  the  net 
movement  of  liquid  water  is  currently  neglected; 

(10)  the  catalyst  layer  is  assumed  to  be  a  thin  interface 
where  sink  and  source  terms  for  the  reactants  and 
enthalpy  are  specified; 

(11)  electroneutrality  prevails  inside  the  membrane.  The 
proton  concentration  in  the  ionomer  is  assumed  to  be 
constant  and  equal  to  the  concentration  of  the  fixed 
sulfonic  acid  groups. 

Most  of  the  above  assumptions  are  a  “standard”  feature 
of  almost  all  previous  modelling  studies.  In  addition,  com¬ 
monly  used  assumptions  in  previous  studies  were:  (i)  one  or 
two-dimensionality;  (ii)  de-coupling  of  ME  A  and  flow 
channel  transport;  (iii)  isothermal  conditions.  In  the  present 
model  these  limiting  assumptions  are  removed.  One  of  the 
other  major  limitations,  namely  the  assumption  of  a  single 
phase  and  of  non-interacting  water  and  gas  phases  in  the 
pores  of  the  GDEs  will  be  addressed  in  future  work  with  the 
implementation  of  a  two-phase  model  with  both  phases 
interacting  in  the  same  computational  domain. 


2.2.  Modelling  domain  and  geometry 

The  computational  domain  that  was  employed  for  the 
simulations  is  shown  in  Fig.  2.  In  addition  to  transport  across 
the  ME  A  (y-direction),  the  formulation  allows  us  to  account 
for  and  investigate  the  effect  of  non-uniform  transfer  rates 
and  species  concentration  along  the  flow  channels  (v-direc- 
tion),  as  well  as  the  three-dimensional  effects  in  the  trans¬ 
verse  ^-direction  due  to  the  geometry  (alternating  open  flow 
channels  with  land  areas).  These  effects  are  expected  to  be 
particularly  important  in  the  regions  of  the  GDEs  under  the 
collector  plates  and  not  directly  exposed  to  the  flow  fields. 
We  take  advantage  of  the  geometric  periodicity  of  the  cell  in 
order  to  reduce  the  size  of  the  computational  domain  and 
hence  computational  cost.  Symmetry  is  therefore  assumed  in 
the  middle  planes  of  the  flow  channel  and  land  areas  and 
hence  only  half  of  each  needs  to  be  incorporated  in  the 
domain.  This  is  a  valid  assumption  as  long  as  no  cross-flow 
takes  place  between  adjacent  channels,  as  in  inter-digitated 
designs,  and  as  long  as  the  region  considered  consists  of 
parallel  straight  channels,  as  is  the  case  for  the  bulk  of  the 
collector  plates  of  an  actual  fuel  cell. 

In  order  to  effectively  implement  the  numerical  solution 
of  the  various  transport  equations,  three  subdomains  were 
defined  within  the  main  computational  domain. 

•  The  main  domain  accounts  for  the  flow,  heat  and  mass 
transfer  of  the  reactant  gases  inside  the  flow  channels  and 
the  GDEs. 

•  Subdomain  I  consists  of  the  MEA  only,  and  accounts  for 
the  heat  flux  through  the  solid  matrix  of  the  GDEs  and  the 
membrane.  Hence,  the  only  variable  of  interest  here  is  the 
temperature.  Exchange  terms  between  this  subdomain 
and  Domain  I  account  for  the  heat  transfer  between  the 
solid  phase  and  the  gas  phase. 

•  Subdomain  II  is  used  to  solve  for  the  flux  of  liquid  water 
through  the  MEA.  The  flux  of  the  water  in  the  membrane 
is  coupled  to  the  electrical  potential  calculated  in  Domain 
III  via  the  Schlogl  equation. 

•  Subdomain  III  consists  of  the  membrane  only  and  is  used 
to  calculate  the  electrical  potential  inside  the  membrane. 


Fig.  2.  Three-dimensional  computational  domain. 
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Fig.  3.  Computational  mesh  for  the  MEA  and  gas-flow  channel  regions. 
The  “empty”  blocks  on  the  left  corresponds  to  the  bipolar  plates,  which 
have  been  left  out  for  clarity. 


The  main  domain  and  three  subdomains  are  all  coupled 
through  appropriate  boundary  exchange  terms  and  an  itera¬ 
tive  solution  procedure.  The  computational  mesh  is  shown  in 
Fig.  3.  Due  to  the  large  number  of  coupled  equations  that  are 
solved  and  the  overall  complexity  of  the  problem,  the 
number  of  computational  cells  was  limited  to  roughly 
80,000.  The  calculations  presented  here  have  all  been 
obtained  on  a  PC  with  a  450  MHz  Pentium  II.  The  number 
of  iterations  required  to  obtain  converged  solutions  ranged 
from  2000  for  lower  current  densities  to  20,000  for  higher 
current  densities;  the  latter  required  about  50  h  of  CPU  time. 

2.3.  Model  equations 


2.3.1.  Notation 

In  the  following,  the  subscript  “g”  denotes  properties  of 
the  gas  phase,  whereas  “1”  stands  for  the  liquid  phase  and 
“s”  for  the  solid  phase.  Different  species  are  denoted  by  the 
subscript  i.e.  the  subscript  “gz”  denotes  the  species 
in  the  gas  phase,  “w”  is  used  for  water  (species),  “sat” 
means  saturation  value.  Cathode  side  and  anode  side  proper¬ 
ties  are  denoted  by  the  subscripts  “c”  and  “a”,  respectively. 


2.3.2.  Fuel  cell  channels 

In  the  fuel  cell  channels,  the  gas-flow  field  is  obtained  by 
solving  the  steady- state  Navier-Stokes  equations,  i.e.  the 
continuity  equation: 

v  •  (pgMg)  =  o  (4) 


and  momentum  equation 

V  •  (PgWg  0  Mg  -  tigX7ug)  =  -V  •(/?  +  §  fig V  •  Mg) 

+  v  •  ^g(Vng)r; 


The  temperature  field  is  obtained  by  solving  the  convective 
energy  equation 

V  •  (pgugHg  -  AgV7g)  =  0  (6) 


here  pa  is  the  gas  phase  density,  u  =  (m,  v,  w)  the  fluid 
velocity  vector,  p  the  pressure,  T  the  temperature,  p  is  the 
molecular  viscosity  and  2  is  the  thermal  conductivity. 


The  total  enthalpy  H  is  determined  from  the  static  (ther¬ 
modynamic)  enthalpy  h  via: 

Hg=  hg+  (7) 

where  the  bulk  enthalpy  is  related  to  the  mass  fraction  y 
and  the  enthalpy  of  each  gas  by: 

hg  =  ygihgi  (8) 

The  mass  fraction  of  the  different  species  obeys  a  transport 
equation  of  the  same  form  as  the  generic  advection-diffusion 
equation: 

V  •  (P&u%y%i)  -  v  •  (/9g£>g,'Vyg()  =  Sgi  (9) 

where  the  i  represents  oxygen  at  the  cathode  side  and 
hydrogen  at  the  anode  side  and  Dgi  is  the  diffusivity  of 
the  species  in  the  background  fluid.  The  source  term  Sgi  is 
determined  by  solving  the  Stefan-Max  well  equations, 
which  account  for  the  diffusion  of  multiple  species  [15]: 


where  v*  is  the  diffusion  velocity  vector  of  species  /,  x  the 
molar  fraction  and  Dy  the  binary  diffusivity  of  any  two 
species. 

The  second  species  on  both  sides  is  water  vapour,  which  is 
assumed  to  exist  at  the  saturation  pressure,  so  that  the  molar 
fraction  of  water  vapour  is  given  by 


(11) 


The  saturation  pressure  of  water  vapour  has  been  approxi¬ 
mated  by  [14] 

log10/4at(r)  =  — 2.1794+0.02953T— 9.1837T2  +  1.4454T3 

(12) 


where  T  is  the  temperature  in  K.  The  mass  fraction  is  related 
to  the  molar  fraction  by 


y& 


Tg jMi 

/A  xglMt 


(13) 


where  M  is  the  molecular  weight  of  the  different  species. 
The  ideal  gas  assumption  leads  to 


(14) 


and  the  bulk  density  becomes: 


1 


Pgi 


The  sum  of  all  mass  fractions  is  equal  to  unity 
=  i 


(15) 


(16) 


which  determines  the  mass  fraction  of  the  third  species  on 
both  sides  (nitrogen  at  the  cathode  and  carbon  dioxide  at  the 
anode). 
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2.3.3.  Gas  diffusion  electrodes 

In  the  porous  GDEs  the  Navier-Stokes  equations  have  to 
be  corrected  for  the  porosity  e  of  the  carbon  fiber  paper. 
Thus,  the  conservation  equation  for  mass  becomes 


v  •  (p  g£g«g)  =  °  (l7> 

whereas  the  momentum  equation  reduces  to  Darcy’s  law: 


Gcyltcy  - 


(18) 


where  kv  denotes  the  hydraulic  permeability.  It  was  men¬ 
tioned  earlier  that  the  liquid  water  pores  are  de-coupled  from 
the  gas  pores,  and  Darcy’s  law  is  again  used  for  liquid  water 
transport: 


(19) 


The  mass  transport  equation  in  porous  media  becomes 


^  ’  (Pg^gWgJg/)  V  •  (PgDgiSgVygi)  Sgi  (20) 

and  the  Stefan-Maxwell  equations  remain  the  same: 


V  *  Xg(  — 


X&Xgj 

Df 


(21) 


where  the  binary  diffusivities  D^f  have  been  corrected  for 
the  porosity.  This  was  done  by  applying  the  so-called 
Bruggemann  correction: 


ne ff  —  n..  *  o1-5 
uij  ~  uu  *  £g 


(22) 


whereas  the  sink  term  for  hydrogen  is  specified  as 


V  =  (26) 

where  F  is  the  Faraday  constant  (96487  C/mol)  and  i  is  the 
local  current  density.  Again,  M  is  the  molecular  weight  of 
the  species  i.  The  product  water  is  assumed  to  be  in  liquid 
form,  and  hence  the  source  term  can  be  written  as 


c  _  ^h2o  . 
6h20(1)  —  2 p  lc 


(27) 


The  generation  of  heat  in  the  cell  is  due  to  entropy 
changes  as  well  as  irreversibilities  due  to  the  activation 
overpotential  f]aci  [11]: 


4c 


T(-ASC) 

nP-F 


+  4 act 


(28) 


where  T  is  the  temperature,  ASC  is  the  entropy  change  in  the 
chemical  reactions,  nQ-  is  the  number  of  electrons  trans¬ 
ferred  and  rjact  is  the  activation  overpotential.  Because  both 
these  contributions  are  small  at  the  anode  side,  the  source 
term  is  currently  neglected  here.  The  overpotential  rjact  can 
be  estimated  a  priori.  This  does  not  introduce  a  large  error, 
since  the  range  of  the  activation  overpotential  at  the  cathode 
side  is  well  known,  depending  on  the  catalyst  loading  and  the 
expected  exchange  current  density. 

As  can  be  seen  from  Eqs.  (25)-(28),  for  an  accurate 
solution  of  the  reactant  gas  distribution  inside  the  fuel  cell, 
it  is  crucial  to  obtain  the  local  current  density  distribution  i, 
which  is  described  by  the  Butler- Volmer  equation  [1]: 


Finally,  the  energy  equation  in  the  diffusion  layer  is  given 
by 

^  '  (pg£gwg^g  —  2g VTg)  =  Sgf{Tg  —  Ts)  (23) 

where  the  term  on  the  right  hand  side  accounts  for  the  heat 
transfer  from  the  solid  matrix  to  the  gas  phase,  f  is  a 
modified  heat  transfer  coefficient  that  accounts  for  the 
convective  heat  transfer  in  W/m2  and  the  specific  surface 
area  m2/m3  of  the  porous  medium.  Hence,  the  unit  of  ft  is 
[W/m3]. 

In  the  solid  matrix  of  the  gas- diffusion  layer,  heat  transfer 
is  calculated  via: 


V  •  (2SVTS)  =  Sgp(Ts  -  Tg) 


(24) 


2.3.4.  Catalyst  layer 

As  mentioned  before,  the  catalyst  layer  is  treated  as  a  thin 
interface,  where  sink  and  source  terms  for  the  reactants  are 
implemented.  Due  to  the  infinitesimal  thickness,  the  source 
terms  are  actually  implemented  in  the  last  grid  cell  of  the 
porous  medium.  At  the  cathode  side,  the  sink  term  for 
oxygen  is  given  by 
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where  c  denotes  the  concentration  of  the  reactants,  and  aa 
and  ac  are  the  so-called  transfer  coefficients.  The  reference 
exchange  current  density  ?Qef  depends  on  various  parameters 
such  as  operating  temperature  and  catalyst  loading,  and  a 
number  of  experiments  have  been  conducted  to  quantify  this 
dependence  empirically  [13,16]. 


2.3.5.  Membrane  model 

In  the  membrane,  the  model  accounts  for  the  flux  of  liquid 
water  and  the  protonic  flux  together  with  the  distribution  of 
electrical  potential  ®  are  being  considered.  The  transport  of 
liquid  water  through  the  membrane  is  governed  by  a  mod¬ 
ified  version  of  Darcy’s  law,  the  Schlogl  equation  [3]: 
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where  k $  denotes  the  hydraulic  permeability,  Zf  is  the  fixed- 
charge  number  in  the  membrane,  Cf  is  the  fixed-charge 
concentration  and  p]  is  the  liquid  water  viscosity.  This 
equation  accounts  for  two  different  water  transport  pro¬ 
cesses:  the  electro-osmotic  drag  whereby  hydrogen  protons 
migrating  through  the  membrane  drag  water  molecules  with 
them,  and  pressure  driven  flux,  which  is  usually  directed 
from  the  cathode  side  to  the  anode  side.  Strictly  speaking,  a 
diffusive  term  has  to  be  accounted  for  as  well,  since  the  back 
diffusion  of  water  can  play  an  important  role  for  water 
management  schemes.  However,  when  the  membrane  is 
fully  humidified,  as  is  assumed  here,  this  term  drops. 

Heat  transfer  in  the  membrane  is  governed  by 

i2 

V2mem  .  vr  =  -  (32) 

K 

where  k  is  the  electric  resistance  of  the  membrane  material. 
Note  that  this  implies  that  the  transport  of  energy  by  the 
liquid  water  is  neglected  and  the  membrane  is  considered  as 
a  heat-conducting  solid.  The  term  on  the  right  hand  side  of 
the  equation  denotes  Ohmic  heating  inside  the  membrane. 
The  local  current  density  inside  the  membrane  is  obtained 
from  Ohm’s  law: 

i  =  -*cV< P  (33) 

Finally,  for  the  electrical  potential  in  the  membrane,  it  can  be 
shown  that  for  a  fully  hydrated  membrane  it  is  governed  by 
the  Laplace  equation  [4]: 

V2<2>  =  0  (34) 


2.3.6.  Bipolar  plates 

Bipolar  plates  serve  to  transfer  the  electrons  and  separate 
different  cells.  Since  the  electrical  conductivity  of  graphite 
2gr  is  high,  Ohmic  losses  are  neglected,  and  the  energy 
equation  reduces  to 

V2gr  •  vr  =  0  (35) 


2.3.7.  Cell  potential 

The  cell  potential  E  is  obtained  by  substracting  all  over¬ 
potentials  (losses)  from  the  equilibrium  potential,  i.e. 

E  ^act  ^mem  (36) 

where  Ejp  is  the  equilibrium  potential  for  a  given  tempera¬ 
ture  and  pressure,  ?7act  the  activation  overpotential  at  both 
sides,  f] Q  are  the  Ohmic  losses  in  the  GDE  plus  contact 
resistances  and  is  the  Ohmic  loss  in  the  membrane. 

The  equilibrium  potential  E®p  is  a  function  of  pressure 
and  temperature  and  is  determined  using  the  Nernst  law  [5], 
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where  p  denotes  the  partial  pressure  of  the  species,  n  is  the 
number  of  electrons  transferred  in  the  reaction,  and  AS  is  the 


change  in  entropy.  The  standard  reversible  potential  E°  is 
given  by  [5] 


o  AG° 

E°  = -  (38) 

nF  v  7 

where  AG°  is  the  change  in  the  standard  Gibbs  free  energy  of 
the  reaction.  Using  standard  values  for  the  entropy  produc¬ 
tion,  Eq.  (37)  yields  [2] 


E°Tp  =  1.229  -  0.83  x  10“3(T  -  298.15)+ 


x  4.31  x  io_5r 


Inpm  +  -  \np0l 


(39) 


2.4.  Boundary  conditions 


Boundary  conditions  have  to  be  applied  for  all  variables  of 
interest  in  each  computational  domain.  In  order  to  reduce 
computational  cost,  we  take  advantage  of  the  geometric 
periodicity  of  the  cell.  Hence  symmetry  is  assumed  in  the  z- 
direction,  i.e.  all  gradients  in  the  z-direction  are  set  to  zero  at 
the  x-y  plane  boundaries  of  the  domain.  With  the  exception 
the  channel  inlets  and  outlets,  a  zero  flux  condition  is  applied 
at  all  x-boundaries  (y-z  planes). 

The  inlet  values  at  the  anode  and  cathode  are  prescribed 
for  the  velocity,  temperature  and  species  concentrations 
(Dirichlet  boundary  conditions).  The  inlet  velocity  is  a 
function  of  the  desired  average  current  density  /ave,  the 
geometrical  area  of  the  membrane  Amea,  the  channel 
cross-section  area  Ac h,  and  the  stoichiometric  flow  ratio  £, 
according  to: 
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where  R  is  the  universal  gas  constant,  7jn  is  the  inlet 
temperature  and  pm  is  the  inlet  pressure. 

At  the  outlets  of  the  gas-flow  channels,  only  the  pressure 
is  being  prescribed  as  the  desired  electrode  pressure;  for  all 
other  variables,  the  gradient  in  the  flow  direction  (x)  is 
assumed  to  be  zero  (Neumann  boundary  conditions). 

Since  the  fluid  channels  are  bordered  by  the  collector 
plates,  no  boundary  conditions  have  to  be  prescribed  here 
and  conjugate  heat  transfer,  impermeability  and  no- slip 
conditions  are  implemented  implicitly  at  solid-fluid  inter¬ 
faces  within  the  domain.  At  the  outer  boundaries  of  the 
bipolar  plates  (y-direction),  boundary  conditions  need  only 
to  be  given  for  the  energy  equation.  This  can  be  done  by 
prescribing  the  heat  flux  or  the  temperature  distribution. 
In  the  present  simulations  a  zero  heat  flux  condition  was 
used,  dT/dy  =  0.  This  is  not  an  entirely  physical  condition, 
but  is  adequate  for  the  present  simulations  in  which  the  focus 
is  on  model  validation  and  identification  of  key  transport 
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processes.  In  future  work  we  intend  to  implement  a  cooling 
flow  channel  in  the  computational  domain  to  remove  arbi¬ 
trary  specification  of  the  heat  flux  or  temperature. 


2.4.1.  Subdomain  I 

As  discussed  earlier,  three  subdomains  are  defined  for 
numerical  efficiency.  Subdomain  I  is  used  to  solve  for  the 
heat  flux  through  the  solid  matrix  of  the  porous  GDE.  The  only 
variable  of  interest  here  is  the  temperature,  and  it  is  assumed 
that  all  the  heat  transfer  takes  place  to  and  from  the  gas  phase 
inside  the  porous  medium  only.  Adiabatic  boundary  condi¬ 
tions  are  therefore  applied  at  all  boundaries  of  this  domain 


At  the  interfaces  exposed  to  the  land  area,  no-flux  conditions 
are  imposed. 


2.4.3.  Subdomain  III 

Finally,  for  the  membrane  domain,  where  the  only  vari¬ 
able  of  interest  is  the  electric  potential,  a  reference  potential 
value  of  zero  is  arbitrarily  set  at  the  anode  side 

0  =  0  (45) 


and  at  the  cathode  side,  the  potential  distribution  at  the 
membrane/catalyst  interface  is  computed  from  [4] 


d0 

dy 


i  F CfVGiinem 

^eff 


(46) 


2.4.2.  Subdomain  II 

For  the  liquid  water  transport  through  the  MEA  in  Sub- 
domain  II,  the  pressure  is  given  at  the  outer  boundaries  of  the 
GDE,  i.e.  the  channel/GDE  interface 

Pa,i  =  3  atm  (43) 

and 

Pci  =  5  atm  (44) 


Table  1 

Geometrical,  operational,  electrode  and  membrane  parameters  for  the  base 
case 


Property 


Value 


Channel  length,  1 
Channel  height,  h 
Channel  width,  wc 
Land  area  width,  w\ 

Electrode  thickness,  te 
Membrane  thickness,  fmem 
Inlet  fuel  and  air  temperature,  T-m 
Air  side  pressure,  pc 
Fuel  side  pressure,  pa 
Air  stoichiometric  flow  ratio,  £c 
Fuel  stoichiometric  flow  ratio,  £a 
Relative  humidity  of  inlet  gases,  2 
Oxygen/nitrogen  ratio,  f 
Gas  phase  electrode  porosity,  eg 
Electronic  conductivity,  o 
Effective  thermal  conductivity,  /cff 
Heat  transfer  coefficient  between  solid 
and  gas  phase,  f 

Transfer  coefficient,  anode  side,  aa 
Transfer  coefficient,  cathode  side,  ac 
Ref.  exchange  current  density,  anode,  i\fd 
Ref.  exchange  current  density,  cathode, 
Entropy  change  of  cathode  side  reaction,  AS?t 
Ionic  conductivity  of  the  membrane,  k 
Protonic  diffusion  coefficient,  Du+ 
Fixed-charge  concentration,  c f 
Fixed- site  charge,  Zf 
Electrokinetic  permeability,  k $ 

Hydraulic  permeability,  kv 

Thermal  conductivity  of  the  membrane,  X 


0.05  m 
l.Ox  10-3  m 
l.Ox  10-3  m 
l.Ox  10-3  m 
0.26x  10-3  m 
0.23  xlO-3  m 
80  °C 
5  atm 
3  atm 
3 
3 

100% 

0.79/0.21 
0.4  [4] 

6000  S/m 
75  W/mK  [9] 

7.0 xlO5  W/m3 

0.5 
1  [IV] 

0.6  A/cm2 

4.4  x  10-7  A/cm2 

-326.36  J/(mol  K)  [11] 

0.06  S/cm 

4.5 xlO-9  m2/s2  [4] 

1200  mole/m3  [4] 

-1 

2.0  x  10_2°  m2  [4] 
1.8xl0-18  m2  [4] 

0.67  W/mK  [10] 


2.5.  Modelling  parameters 


The  parameters  used  for  the  base  case  simulations  pre¬ 
sented  here  are  shown  in  Table  1 .  Since  the  model  presented 
here  was  initially  developed  by  extending  the  formulation  of 
Bernardi  and  Verbrugge  [3,4]  to  three-dimensions,  many  of 
the  key  parameters  defined  by  these  authors  are  still  used 
here.  It  is  important  to  note  that  because  this  model  accounts 
for  all  major  transport  processes  and  the  modelling  domain 
comprises  all  the  elements  of  a  complete  cell,  no  parameters 
needed  to  be  adjusted  in  order  to  obtain  physical  results. 
When  comparing  our  results  with  experiments  published  in 
the  literature,  however,  many  of  the  experimental  data  such 
as  the  stoichiometric  flow  ratio  or  the  exact  cell  geometry 
and  dimensions  are  unknown,  which  makes  a  quantitative 
comparison  difficult.  The  strength  of  the  current  model  is 
clearly  to  perform  parametric  studies  and  explore  the  impact 
of  various  parameters  on  the  transport  mechanisms  and  on 
fuel  cell  performance. 

For  the  binary  diffusivities  Dy  required  in  the  Stefan- 
Maxwell  equations,  experimentally  obtained  values  at  atmo¬ 
spheric  pressure  po  were  taken  and  scaled  with  the  tem¬ 
perature  and  pressure  according  to  [6] 


Dij  —  Dij(To,po) 


P_ 

Po 


(47) 


Table  2  lists  the  reference  binary  diffusivities. 


Table  2 


Binary  diffusivities  at  reference  temperatures 


Gas-pair 

Reference 
temperature,  T{)  (K) 

Binary  diffusivity, 
Dtj  (cm2/s) 

^h2-h2o 

307.1 

0.915 

Dr2-co2 

298.0 

0.646 

Dr2o-co2 

307.5 

0.202 

Do2-r2o 

308.1 

0.282 

Do2-n2 

293.2 

0.220 

Dr2  o-n2 

307.5 

0.256 
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3.  Results  and  discussion 

3.1.  Fuel  cell  performance 

The  fuel  cell  performance  is  shown  in  terms  of  the 
polarization  in  Fig.  4.  The  results  agree  well  with  the 
experiments  [17]  in  the  low  and  intermediate  current  density 
region.  The  discrepancy  at  high  current  densities  is  attrib¬ 
uted  to  the  fact  that  experimental  data  at  the  higher  current 
densities  was  insufficient,  and  in  fact  the  experimental  curve 
in  this  region  is  a  curve  fit  weighted  by  the  lower  current 
density  data.  It  should  also  be  pointed  out  that  the  exact 
geometry  of  the  fuel  cell  used  in  the  experiments  is  not 
known.  The  ability  of  the  present  model  to  reproduce  the 
polarization  curve  is  a  necessary  validation  check  but  by  no 
means  especially  informative  since  one  can  always  obtain 
good  agreement  between  experimental  results  and  any 
model  that  somehow  captures  the  logarithmic  drop-off  in 
the  low  current  density  region  and  the  Ohmic  losses  inside 
the  membrane.  The  strength  of  the  numerical  approach  is  in 
providing  detailed  insight  into  the  various  transport  mechan¬ 
isms  and  their  interaction,  and  in  the  possibility  of  perform¬ 
ing  parameters  sensitivity  analyses. 


0.0035 


0.05 


0.0020  0.00 


XH2 

0.58 

0.56 

0.54 

0.52 

0.50 

0.48 

0.46 

0.44 


X02 

020 

0.18 

0.16 

0.14 

0.12 

0.10 

0.08 

0.06 


Fig.  5.  Three-dimensional  distribution  of  reactant  gases  in  the  gas-flow 


3.2.  Reactant  gas  and  temperature  distribution 

One  of  the  major  advantages  of  using  a  model  as  detailed 
as  the  one  presented  here  is  the  detailed  distribution  of  the 
reactant  gases  inside  the  fuel  cell.  Such  distributions,  which 
cannot  be  measured  in  situ,  provide  valuable  information 
about  the  onset  of  concentration  losses  and  their  effect  on  the 
limiting  current  density.  Fig.  5  shows  the  reactant  gas 
distribution  in  the  fuel  cell  channels  and  the  gas-diffusion 
layers.  The  molar  fraction  of  oxygen  decreases  noticeable 
inside  the  gas-diffusion  layer,  and  the  effect  of  oxygen 
depletion  is  significant,  particularly  under  the  land  areas. 
A  realistic  prediction  of  the  limiting  current  density  cannot 
be  obtained  from  a  two-dimensional  model  in  which  infi¬ 
nitely  wide  channels  are  assumed  and  the  area  that  is  not 
exposed  to  the  gas-flow  channels  is  not  accounted  for. 


Fig.  4.  Comparison  of  the  experimental  and  simulated  polarization  curves. 
Also  shown  is  the  power  density  curve  obtained  with  the  model. 


The  effect  of  the  “land  area”  is  even  more  pronounced  at 
higher  current  densities,  as  shown  in  Fig.  6.  Furthermore,  for 
a  higher  current  density  at  otherwise  similar  conditions,  the 
oxygen  molar  fraction  at  the  catalyst  layer  decreases  due  to 
diffusion  limitations.  Note  that  the  stoichiometric  flow  ratio 
is  the  same  for  both  cases,  i.e.  three  times  the  amount  of 
oxygen  consumed  enters  the  cell.  Since  the  composition  of 


Fig.  6.  Molar  oxygen  fraction  at  the  catalyst  layer  for  two  different  current 
densities:  0.2  A/cm2  (upper)  and  1.0  A/cm2  (lower). 
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Fig.  7.  Dimensionless  current  density  distribution  at  the  catalyst  layer 
interface  for  three  different  nominal  current  densities:  0.2  A/cm2  (upper); 
0.8  A/cm2  (middle);  1.4  A/cm2  (lower).  The  dashed  lines  represent  the 
boundaries  between  the  flow  channel  and  the  land  area. 

the  inlet  gases  is  the  same,  this  is  accomplished  by  an 
increase  in  velocity  for  a  higher  current  density.  At  an 
average  current  density  of  1.4  A/cm2,  the  molar  oxygen 
fraction  at  the  catalyst  layer  is  almost  zero  throughout  the 
interface,  indicating  that  the  limiting  current  density  has 
been  reached. 

The  detailed  oxygen  distribution  at  the  catalyst  layer  is  of 
importance,  because  it  determines  the  local  current  density, 
according  to  Eq.  (29).  Assuming  a  constant  activation  over¬ 
potential  throughout  the  catalyst  interface,  this  results  in  local 
current  densities  as  shown  in  Fig.  7.  The  local  current 
densities  here  have  been  normalized  by  the  nominal  current 
density  in  each  case.  It  can  be  seen  that  for  a  low  nominal 
current  density,  the  distribution  is  quite  uniform  with  varia¬ 
tions  from  75  to  125%  relative  to  the  nominal  current  density 
of  0.2  A/cm2 .  This  changes  for  intermediate  current  densities, 
where  under  the  land  areas  a  noticeable  decrease  takes  place 
and  the  minimum  current  fraction  drops  below  40%  of  the 
average  (nominal)  current  density  of  0.8  A/cm2.  The  max¬ 
imum,  however,  remains  almost  the  same  at  just  over  130%. 
This  pattern  changes  further  at  high  current  densities,  where 
the  maximum  local  current  density  can  be  as  high  as  three 
times  the  average  current  density  near  the  cathode  side  inlet. 

With  an  increase  in  current  density,  the  fraction  of  the 
current  density  generated  under  the  area  exposed  to  the 


Fig.  8.  Fraction  of  the  total  current  generated  under  the  area  exposed  to  the 
gas-flow  channels. 

gas-flow  channels  increases  steadily,  as  shown  in  Fig.  8.  The 
fraction  of  the  current  generated  under  the  channel  area 
increases  almost  linearly  with  the  current  density,  the  max¬ 
imum  being  just  under  80%  at  a  current  density  close  to  the 
limiting  current  density.  This  is  important,  since  it  provides 
an  indication  of  how  well  the  catalyst  is  being  used  in  these 
areas.  For  optimal  fuel  cell  performance,  a  uniform  current 
density  generation  is  desirable,  and  this  could  only  be 
achieved  with  a  non-uniform  catalyst  distribution,  possibly 
in  conjunction  with  non-homogeneous  GDEs. 

One  of  the  most  important  features  distinguishing  the 
present  model  is  the  fact  that  it  is  non-isothermal.  The 
temperature  distribution  inside  the  fuel  cell  has  important 
effects  on  nearly  all  transport  phenomena,  and  knowledge  of 
the  magnitude  of  temperature  increases  due  to  irreversibil¬ 
ities  might  help  preventing  failure.  Fig.  9  shows  that  the 
increase  in  temperature  can  exceed  several  degrees  Kelvin 
near  the  inlet  area,  where  the  local  current  density  is  highest. 
Due  to  the  low  electric  conductivity  of  the  polymer  electro¬ 
lyte,  the  temperature  maximum  occurs  inside  the  membrane. 
In  general,  the  temperature  at  the  cathode  side  is  slightly 
higher  than  at  the  anode  side;  this  is  due  to  the  reversible  and 
irreversible  entropy  production.  The  detailed  temperature 
distribution  is  also  key  for  the  eventual  extension  of  the 
current  model  to  include  multiphase  phenomena.  Phase 
change  is  not  yet  accounted  for,  but  the  mechanism  shall 
be  briefly  outlined  here.  In  the  first  place,  the  temperature 
rise  at  the  cathode  determines,  how  strong  the  under- satura¬ 
tion  of  the  gas  phase  in  this  area  is.  This  in  turn  gives  rise  to 
the  evaporation  of  liquid  product  water,  which  provides 
cooling  and  thus  offsets  the  temperature  rise.  Hence,  any 
implementation  of  phase  change  in  an  isothermal  model  can 
not  lead  to  physically  representative  results. 

Finally,  the  water  flux  and  the  potential  distribution  is 
presented  in  Fig.  10.  The  flux  of  liquid  water  is  governed  by 
three  mechanisms:  electro-osmotic  drag  from  the  anode  to 
the  cathode,  back  diffusion  driven  by  a  concentration  gra¬ 
dient  from  the  cathode  to  the  anode  and,  if  a  pressure 
gradient  exists,  convection  which  is  usually  directed  from 
the  cathode  to  the  anode  in  order  to  counter-balance  the 
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Fig.  9.  Three-dimensional  temperature  distribution  inside  the  fuel  cell  at  a 
nominal  current  density  of  1.2  A/cm2. 


Fig.  11.  Net  drag  coefficient  a  for  the  flux  of  water  through  the  membrane 
for  two  different  values  of  the  electrokinetic  permeability  of  the 
membrane. 


electro-osmotic  drag.  Assuming  a  fully  humidified  mem¬ 
brane  implies  no  concentration  gradient  exists,  and  in  any 
case,  the  effect  of  the  diffusion  on  the  water  flux  through  the 
membrane  has  been  found  to  be  small  compared  to  convec¬ 
tion  and  electro-osmotic  drag.  As  illustrated  in  Fig.  10,  the 
pressure  gradient  might  outweigh  the  effect  of  the  electro- 
osmotic  drag  for  low  current  densities,  and  the  net  flux  of 
water  is  directed  from  the  cathode  to  the  anode.  Humidifica¬ 
tion  schemes  have  to  be  therefore  devised  for  the  cathode 
side  only.  At  a  current  density  of  0.4  A/cm2,  the  flux  of  water 
is  reversed.  Near  the  inlet  area,  where  the  local  current  is 


Fig.  10.  Liquid  water  flux  (vectors)  and  potential  distribution  (contours) 
inside  the  membrane  for  three  different  current  densities:  0.2  A/cm2 
(upper);  0.4  A/cm2  (middle);  0.6  A/cm2  (lower). 


strongest,  the  electro-osmotic  drag  dominates  convection, 
and  water  flows  from  the  anode  to  the  cathode.  Near  the 
outlet,  the  current  is  weaker  and  the  net  water  flux  is  directed 
towards  the  anode. 

The  net  flux  of  water  through  the  membrane  is  often 
characterized  by  the  parameter  a,  the  net  amount  of  water 
molecules  dragged  through  the  membrane  per  hydrogen 
proton.  Fig.  1 1  shows  the  net  water  flux  for  two  different 
values  of  the  electrokinetic  permeability  of  the  membrane.  If 
the  value  cited  by  Bernardi  and  Verbrugge  is  used  [4],  a 
becomes  unphysically  high.  A  comparison  with  the  model 
published  by  Yi  and  Nguyen  [21]  reveals  that  reducing  the 
permeability  to  2.0  x  10  20  m2  as  was  used  in  the  current 
base  case  yields  a  more  realistic  values  for  a.  Still,  these 
results  differ  substantially  from  the  experimentally  deter¬ 
mined  values  of  a,  which  range  from  0.6  at  low  current 
densities  to  around  0.3  for  intermediate  current  densities  in 
the  absence  of  a  pressure  gradient.  The  membrane  model  has 
therefore  to  be  improved  in  order  to  predict  the  amount  of 
water  that  need  to  be  supplied  at  the  electrodes  in  order  to 
prevent  drying  out  of  the  membrane. 

4.  Conclusions 

A  comprehensive  three-dimensional  computational 
model  of  a  PEM  fuel  cell  has  been  developed.  With  the 
exception  of  phase  change,  the  model  accounts  for  all  major 
transport  phenomena  in  the  flow  channels,  electrodes  and 
electrolyte  membrane.  Results  that  are  physically  consistent 
and  in  good  agreement  with  available  experimental  data  are 
obtained.  The  three-dimensional  nature  of  the  distribution  of 
flow  velocities,  species  concentration,  mass  transfer  rates, 
electric  current  and  temperature  was  clearly  illustrated  by 
the  simulations.  The  capabilities  of  the  model  for  providing 
detailed  insight  into  water  transport  mechanisms  and  the 
onset  of  mass  transport  limitations  was  demonstrated,  and  its 
potential  for  parametric  studies  of  interest  in  design  and 
prototyping  were  also  illustrated. 
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Though  this  comprehensive  model  represents  to  our 
knowledge  the  only  three-dimensional  model  developed 
to  date,  a  number  of  important  mechanisms  and  phenomena 
need  to  be  incorporated  to  increase  its  generality  and 
usefulness.  Work  is  currently  in  progress  to  account  for 
(a)  the  effect  of  phase  change  of  water,  and  (b)  partial 
dehydration  of  the  membrane  during  operation.  Finally 
the  electrochemistry  in  the  computational  model  relies 
on  first  order  kinetics  and  empirical  data.  Further  work  in 
this  area  is  required  on  both  theoretical  and  experimental 
fronts. 
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